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The Functional Measurement approach, proposed within the theoretical 
framework of Information Integration Theory (Anderson, 1981, 1982), can be 
a useful multi-attribute analysis tool. Compared to the majority of statistical 
models, the averaging model can account for interaction effects without 
adding complexity. The R-Average method (Vidotto & Vicentini, 2007) can 
be used to estimate the parameters of these models. By the use of multiple 
information criteria in the model selection procedure, R-Average allows for 
the identification of the best subset of parameters that account for the data. 
After a review of the general method, we present an implementation of the 
procedure in the framework of R-project, followed by some experiments 
using a Monte Carlo method. 


Multi-attribute models generally follow three steps: evaluation of the 
attributes, integration of the obtained subjective dimensions, and a 
conclusive stage. In the last stage, the results of the previous processes are 
transformed into a ranking order, a set of pairwise preferences or a rating 
over some real interval (Lynch, 1985; Oral & Kettani, 1989). A subset of 
these models, proposed by Anderson (1981, 1982), identifies the averaging 
process as one of the widely used cognitive integration rules. The averaging 
process uses a weight and scale value parameters representation. Ratio 
scales are involved in the measurement of weights and equal-interval scales 
are used to measure values (Zalinski & Anderson, 1989). 

Furthermore, the method of sub-designs, proposed by Norman (1976) 
and Anderson (1982), allows for complete identifiability of these parameters 
by adjoining selected sub-designs to the full factorial design. For instance, a 
full three-way design (A x B x C) can be supplemented with three two-way 
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sub-designs (A x B, A x C and B x C), and with three one-way sub-designs 
(A, B and C). 

In the current literature, this method is applied in different fields of 
Psychology, especially in Social (Falconi & Mullet, 2003; Girard, Mullet, & 
Callahan, 2002; Wang & Yang, 1998), Cognitive (Oliveira et ah, 2006) and 
Developmental (Jager & Wilkening, 2001) Psychology. A computer tool to 
cope with factorial design of Functional Measurement has recently been 
developed by Mairesse, Hofmans, and Theuns (2008). 

Despite the widespread use of the methodology, at present there are 
few tools for the estimation of averaging parameters. The first estimation 
procedure was developed by Zalinski (1984, 1986, 1987) and implemented 
in FORTRAN language. It had the capability to deal with the EAM (Equal 
Averaging Models) and the Complete DAM (Complete Differential 
Averaging Model) and to analyze the data from one subject and one session 
at a time. In addition, it did not account for information criteria (Wang & 
Yang, 1998). 

A second procedure has more recently been developed by Vidotto & 
Vicentini (2007) and focuses on the use of a different minimization 
algorithm for bounded parameters and on the introduction of an information 
criterion able to test several DAMs including the Complete DAM. The 
method has been implemented in the R-project environment: a modern, 
open-source and widely used framework for statistical analysis (R 
Development Core Team, 2009). 


AVERAGING MODELS 

Averaging models involve a two parameters representation: the scale 
value sjh which represents the psychological scale value of the y'-th level of 
an attribute k of some overt dimension, and the weight wjh which represents 
its importance in the integrated response (Anderson, 1965). 

The averaging model represents the integrated response, R, as: 

R = ll^ s i with ° = Z]Io W < (1) 

where the multi-index i = (j,k) accounts for the N overall number of 
attributes that are included in each stimuli set. Also, the condition i = 0 
accounts for prior memorial information (it is also named the prior belief or 
initial state). Note that j = 1, ... J/ c indexes the number of levels of a factor 
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whereas k = 1, ... K indexes the factors of an experimental design. The 
relative weights of each stimulus thus depend on the other stimuli in the set. 

A factor k is said to be equally-weighted if wjk = Wk for every j. Hence, 
the denominator of equation (1) has the same value in each cell of the 
design and can be absorbed into an arbitrary scale unit. If all the factors are 
equally-weighted, then the model is called an Equal-weight Averaging 
Model. 

However, accounting for crossover effects needs the introduction of 
different weight parameters that make the model more complex. These 
families of Differential-weight Averaging Models allow each stimulus (or 
group of stimuli) to have its own weight as well as its own scale value. The 
sum of the absolute weights in the denominator of equation (1) varies 
therefore from cell to cell in the design and the model becomes inherently 
non-linear. Nevertheless, this non-linearity generally introduces some 
analytical and statistical problems with regards to uniqueness, bias, 
convergence, reliability, and goodness of fit (Zalinski & Anderson, 1991). 


R-AVERAGE - THE METHOD 

The R-Average method chooses the optimal model according to the 
“Ockham’s razor”: i.e. the one that fits empirical data by using the smallest 
set of weight parameters. Three goodness-of-fit indexes, based on Residual 
Sum of Squares (RSS) can be used to identify such a model: the adjusted R- 
square, the Akaike Information Criterion (AIC; Akaike, 1974, 1976), and 
the Bayesian Information Criterion (BIC; Schwarz, 1978; Raftery, 1995). 

Starting from the EAM as a baseline, one single weight parameter (or 
a set of them in the following steps) is changed and accepted (or rejected) 
for successive iteration on the result of a comparison between the baseline 
and the new model goodness of fit indexes (whenever ABIC < 2, difference 
in the AIC index is considered). The procedure iterates until no further 
improvements appear and represents a compromise between efficiency and 
performance. It has the capability of providing reliable estimation for each 
trial and for all the repeated measurements on a subject and can account for 
many repetitions to achieve reliable parameter estimation. 

The weight and scale value parameters are estimated by minimizing 
the RSS of the model. This is perfonned with the L-BFGS-B algorithm 
implemented by Byrd et al. (1995), which is an extension of the limited 
memory algorithm L-BFGS. Unlike the latter, the L-BFGS-B method is 
useful for solving large non-linear optimization problems with simple 
bounds on the variables. This algorithm does not require computation of 
second derivatives and the knowledge of the objective function structure. 
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The search direction employs a two-stage approach: the first stage identifies 
a set of active variables using the gradient projection method; in the second 
stage, a quadratic model is approximately minimized with respect to these 
free variables. Once a search direction is established, a line search is 
performed using a method described by More and Thuente (1990). If 
bounds are active, the algorithm stops when the nonn of the projected 
gradient is sufficiently small. 

The choice of a method for solving constrained optimization problems 
stems on the Zalinski's (1987) recommendation that reliable estimations of 
weights can be provided when the minimization function is bounded, but 
other solutions, like Simulated Annealing methods, are also possible. 


R-AVERAGE - THE IMPLEMENTATION 

The R-Average procedure has been implemented as a computer- 
library within the R-project (R Development Core Team, 2009). The R- 
Average library is well integrated in the R framework and it is specifically 
designed to easily manage data with several subjects and repeated 
experiments. The package, which includes R help pages and sample inputs, 
is available from the authors of this paper. 

For command-line beginners a practical graphical interface has been 
developed. 

The algorithm selects the combination of weights that shows the best 
fit indexes by considering each possible DAMs. It can be used for 
estimating weight and scale value parameters of a factorial design both for 
single subject, with repetitions, and for an entire sample; and summarize 
them into tables. It works on any number of factors and levels, and provides 
several ways to analyze data by handling the principal attributes of the 
functions that are listed below: 

• Data: a matrix object containing the experimental data. The first 
column is filled with the initial state values whereas the others 
contain the sub-design response values (in order: one-way sub¬ 
designs, two-way sub-designs, etc.) and the full factorial design. If 
any sub-design is not available, the corresponding matrix columns 
must be loaded with NA. 

• Lev. a vector containing the number of the levels for each factor. 

• Range : a vector containing the range of the scale responses R. 

• Start: a vector containing the scale and the weight values indicating 
the starting point of the computational algorithm. 
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• Lower: a vector containing the lower bound of the scale and weight 
values. 

• Upper: a vector containing the upper bound of the scale and weight 
values. 

• All: a logical attribute that allows the procedure to test all the 
possible subset of DAMs with different weight parameters, or to 
restrict the analysis to a selected subset. 

• Equal, weights: a numeric attribute that allows to fix the number of 
possible equal weights. 

• Delta.weights: a numeric attribute that allows to choose a cut-off 
value at which different weights are considered equal. 

• IC.diff: a vector containing the cut-off values for BIC and AIC at 
which different models with similar BIC/AIC are considered 
equivalent. 

• Verbose: a logical attribute that allows to print general infonnation 
on every step of the information criterion procedure. 

• Maxit: the maximum number of minimization algorithm iterations, 
55 for default. 

• Method: the method followed by the algorithm for estimating the 
parameters. L-BFGS-B is the default option. However, due to the 
fact that results of R-Average implementation sometimes appear to 
be stuck in local minima (since they vary from the ideal values) a 
simulated annealing algorithm (option SANN) can be chosen. 
Nevertheless, its native implementation in the R framework (when 
called by the R-Average routine) has revealed to be time-consuming 
and require further analyses and adaptation. At present we highly 
discourage this option to be run. 


MONTE CARLO SIMULATIONS 

EFFECTS OF PARAMETERS ORDER 

Starting from a strictly increasing monotonic factor A (with levels si < 
S 2 < S 3 ), five simulations were run using different permutations of A for the 
levels of factor B. The schema of the random synthetic data for each of the 
1000 iterations is reported in Table 1. 




Estimated Responses 
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Table 1. Schema of the random synthetic data for each iteration. 


Trend 

1 

Factor A 


Factor B 

1 

Sl 

S2 

S3 

Sl 

S 2 

S3 

2 

Sl 

S2 

S3 

S3 

s 2 

Sl 

3 

Sl 

S2 

S3 

s 2 

Sl 

S3 

4 

S2 

Sl 

S3 

s 2 

Sl 

S3 

5 

S2 

Sl 

S3 

Si 

s 2 

S3 


In every iteration a 5-row data matrix of random values was 
generated, followed by the addition of a normal N(0,1.5) distributed error, in 
order to calculate the responses Rjk. Results show that there are no effects in 
parameters estimation due to the factor’s position in the algorithmic 
procedure or due to the monotonicity of the factors ( R 2 = 0.9778 for all 
trends). The scatter plot is reported in the left panel of Figure 1. 


Theoretical R vs. Estimated R 


Scale value parameters 




Theoretical Responses 


Theoretical s 


Figure 1. On the left: Estimated parameters versus theoretical 
parameters for Trend 1 to 5. On the right: Effects of the bounding on 
extreme value parameters (white: s 2 , light gray: s 2 , dark gray: S 3 ). 

RELIABILITY OF THE PROCEDURE 

Synthetic responses Rjk for factorial design were calculated for a 3 x 
3 experimental design. In the setting of Monte Carlo simulation were 
included different rows number (5,10,15) and different normal error 
(N(0,.5), N(0,1), N(0,1.5)) of the data matrix. These two attributes have 
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been manipulated in order to simulate both different sets of measure 
repetitions and the influence of response variability on the data. 

In Table 3 the results of these Monte Carlo simulations are shown. 
The higher the number of rows, or the lower the standard deviation of the 
error, the better the estimation of the parameters is. 

In particular the procedure has good performance in recognizing the 
appropriate order of the parameters, as can be seen in Table 2. This fact is 
particularly important; in spite of the high variability that sometime affects 
the parameters estimation, at least the necessary condition of ordering is 
often fulfilled, ensuring a solid base for further improvement of the 
estimated values by means of iterative methods, R-Average attributes 
handling or other statistical methods. 


Table 2. Number of correctly estimated parameter order in the MC- 
simulations. 


Rows 

Error 

W >A1< 

WA2 

WA2< 

Wa3 

WBl > 

Wr2 

W/}2 < 

WB3 

Scale 

A 

Scale 

B 



1000 

1000 

1000 

1000 

1000 

1000 

5 

1.0 

985 

993 

954 

909 

979 

863 


1.5 

873 

942 

835 

793 

820 

637 


0.5 

1000 

1000 

1000 

1000 

1000 

1000 

10 

1.0 

1000 

999 

994 

983 

999 

977 


1.5 

980 

986 

955 

890 

966 

845 


0.5 

1000 

1000 

1000 

1000 

1000 

1000 

15 

1.0 

1000 

1000 

999 

996 

1000 

995 


1.5 

995 

996 

982 

937 

991 

919 


Results, however, deserve different considerations depending on the 
type of estimated parameters: 5-type parameters follow a normal distribution 
centered on the true value of the parameter and with a standard deviation 
that increases with the error standard deviation, but decreases with the larger 
sets of available data. Also, 5-parameters that are closer to the bounding 
value (0-20) are generally estimated with an asymmetric distribution as can 
be inferred by the right panel of Figure 1. 

In contrast, vv-type parameters do not follow a nonnal distribution. In 
Table 3 the ratios of the u-paramctcrs are reported, since their proportions 
are more important then their values. 






















Table 3. Results of Monte Carlo simulations. Estimated parameters are reported for each condition on 
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It can be seen that exact estimation of the real value of the w-ratios are less 
precise than those of ^-parameters; in particular, a higher error standard 
deviation increases the difficulty in estimating the exact value. Nevertheless, 
as we noticed before, more reliable values of the ratios can always be 
searched if the fundamental requirement of ordering is fulfilled. 

EFFECTS OF ATTRIBUTES HANDEING 

Monte Carlo simulations have been run to test the procedure 
improvement in parameters estimation with different settings of the 
algorithmic attributes. Five different attributes have been handled in the 
simulation for the conditions of “5 rows, 1.5 variance” and “10 rows, 1.0 
variance”; the first condition was chosen to reproduce an experiment with a 
low number of trials and an high variance in data, the second condition to 
reproduce an experimental situation with a good control in data variance and 
an adequate number of trials. Results are described below and summarized 
in Tables 4, 5 and 6. 

• Changes of bounds or lower and upper attributes: results show that 
changing the range of the weights parameters does not affect the 
reliability of the procedure since the algorithm is based on the 
conservation of the ratios between vv-paramctcrs. The scale for w- 
parameters can thus be chosen according to the experimental setting. 

• Attribute equal.weights set equal to 2: results show that equating a 
priori a couple of parameters could improve estimation reliability 
(see in Tables 4, 5 and 6). 

• Attribute delta.weights set equal to 0.5: results show that setting the 
criterion value to 0.5 at which two weights are considered equal 
makes order recognition more stable (as can be seen in Table 4) but 
worsen the mean estimated parameters value in presence of high 
data variability (5-1.5 simulation, Table 5); instead, in the presence 
of contained variability the results are improved (10-1.0 simulation, 
Table 6). 

• Attributes IC.diff set to 1: using a restrictive acceptance criterion for 
BIC selection slightly worsens the reliability of the procedure in the 
presence of larger variance. Reliability of the order of the parameters 
of factor B is improved, yet the reliability of the scale order of the 
other parameters is reduced. Notice that the change in IC.diff does 
not affect mean estimated parameters as can be seen in Tables 5 and 
6 . 
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• Attributes All set to True: allowing the algorithm to span all the 
possible DAMs can increase the reliability of the procedure. Results 
show that there is a strong improvement both in the exact estimation 
of parameters (Tables 5 and 6) and in the recognition of the correct 
parameters order, as can be seen in Table 4. 

• Combination of All equal to True and delta.weights equal to 0.5: 
performance can than be strongly improved by combining the best 
previous attribute manipulations as in Table 4, 5 and 6. 


Table 4: Results of Monte Carlo simulations for different attributes 


manipulation conditions. 


5 Rows 1.5 Error 

w A i< 

W A 2 

W A 2< 

W A 3 

w B i> 

H>B2 

WB2 < 

W B 3 

Scale A 

Scale B 

Default 

873 

942 

835 

793 

820 

637 

Equal.weights=2 

870 

943 

823 

790 

819 

626 

Delta.weights=0.5 

963 

956 

928 

834 

919 

765 

IC.diff=l 

825 

943 

790 

793 

774 

597 

All=True 

949 

952 

908 

877 

901 

785 

A11=TRUE and 
delta.weights=0.5 

958 

950 

924 

872 

910 

796 

10 rows 1.0 error 

w A i< 

W A 2 

W A 2< 

W A 3 

w BI > 

H>B2 

WB2 < 

WB3 

Scale A 

Scale B 

Default 

1000 

999 

994 

983 

999 

977 

Equal.weights=2 

1000 

999 

993 

985 

999 

978 

Delta.weights=0.5 

1000 

999 

997 

988 

999 

985 

IC.diff=l 

998 

999 

984 

979 

997 

963 

All=True 

1000 

999 

994 

983 

999 

977 

A11=TRUE and 
delta.weights=0.5 

1000 

999 

995 

988 

999 

983 
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GENERAL DISCUSSION 

R-Average consists of a general procedure and an R-library for 
parameter estimation and model selection of the averaging model. It is well 
integrated in the R framework (R Development Core Team, 2009). 

The procedure, based on the “Ockham’s Razor” and on the bounded 
minimization algorithm L-BFGS-B, provides good estimates of the 
parameters of the averaging models and several goodness-of-fit indexes 
useful for model comparison. It also allows the management of several 
repetitions both for single subject and for group data and to estimate 
parameters for incomplete factorial designs. 

In addition, it provides several attributes useful to improve analysis: 
for instance, the number of equal weights can be established a priori and 
several criteria can be set to adapt the procedure flow to the experimental 
needs. These attributes are listed previously in this text. 

Results of Monte Carlo simulations show good reliability of the 
procedure in parameters estimation and excellent reliability in recognition of 
the parameters importance order. This latter result ensures good behavior of 
the procedure and implies that more accurate estimations of the parameters 
could be obtained. For this purpose, further developments of the procedure 
will allow the selection and fixing of the parameter values. Then a Monte 
Carlo simulation study can be conducted based on experimental data, in 
order to obtain more reliable parameters estimation. 

Finally, since reliability reduces with the increase of data variability, 
one remaining issue concerns the theoretical validity of group analyses: this 
issue was previously raised by Zalinski and Anderson (1991) because these 
analyses can cause bias in the estimated weights and scale values; 
nevertheless, if the variability of the responses is not extreme, the R- 
Average procedure can be carefully used even for group data. 
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